﻿// mod09ga_prepare.cpp : 此文件包含 "main" 函数。程序执行将在此处开始并结束。
//使用C盘 OSGeo4W的GDAL库快速开发GDAL程序。
// 最后结论，使用COG，第0级，第1级的问题，不想再弄了就这样吧。采用单瓦片VRT的形式读取数据。 2024/5/15 13:52

#include <string>
#include <vector>
#include <iostream>
#include "gdal_priv.h"
#include "gdalwarper.h"
#include "vrtdataset.h"
#include "../../sharedcodes/wdatetime.h"
using namespace std;


int tryWarpedToMem();
int tryOutputTifFile();
int tryWaredVrt();
int64_t tryATile(string fname, int tilez, int tiley, int tilex);//构建瓦片VRT直接读
int64_t tryATile2(string fname, int tilez, int tiley, int tilex);//构建全局VRT然后读取瓦片
void tryOneFileTiles(string fname);
void tryOneFileTiles2(string fname);

int main(int argc,char* argv[])
{
    cout << "version v1 warp output a tif is ok" << endl;
    cout << "version v2 try warped to a mem dataset. OK." << endl;
    cout << "version v3 try vrt, this worked! GDALDestroyGenImgProjTransformer exception sovleved" << endl;
    cout << "version v4.1 try performance for no-ov-tiff, ov512-tiff, cog-tiff. top-level tiles are too slow about 100ms per tile in 0-2 level." << endl;

    //tryWarpedToMem();
    //tryOutputTifFile();
    //tryWaredVrt();

    string noOv = "E:/data/Try_COG_TilerPerformance/ndvi9_rgb_rendered_no_overview.tif";
    string lineOv = "E:/data/Try_COG_TilerPerformance/ndvi_rgb_rendered_overviewed.tif";
    string cog = "E:/data/Try_COG_TilerPerformance/ndvi9_rgb_rendered_cog.tif";
    string normalTifWithTiled = "E:/data/Try_COG_TilerPerformance/tiled_normal_tif.tif";
    string tiledOvr = "E:/data/Try_COG_TilerPerformance/tiled_ovr_ntif.tif";
    tryOneFileTiles(tiledOvr);
    tryOneFileTiles2(tiledOvr);
    tryOneFileTiles(normalTifWithTiled);
    tryOneFileTiles2(normalTifWithTiled);
    tryOneFileTiles(noOv);
    tryOneFileTiles2(noOv);
    tryOneFileTiles(lineOv);
    tryOneFileTiles2(lineOv);
    tryOneFileTiles(cog);
    tryOneFileTiles2(cog);

//ndvi_rgb_rendered_overviewed.tif
//ndvi9_rgb_rendered_cog.tif
//ndvi9_rgb_rendered_no_overview.tif
    return 0;
}
//return current milliseconds
int64_t currentMilliSeconds()
{
    unsigned long now = std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::system_clock::now().time_since_epoch()).count();
    return now;
}

int64_t tryATile2(string fname, int tilez, int tiley, int tilex)
{
    GDALAllRegister();
    int numXTiles = pow(2, tilez);
    double reso = 360.0 / numXTiles / 256.0;
    double left = -180.0 + tilex * 256 * reso;
    double top = 90.0 - tiley * 256 * reso;

    int64_t ms0 = currentMilliSeconds();
    GDALDataset* inds = (GDALDataset*)GDALOpen(fname.c_str(), GA_ReadOnly);
    int outxsize = 256* numXTiles;
    int outysize = 256* numXTiles/2;
    double trans[6] = { -180.0 , reso , 0 , 90.0 , 0, -reso };
    string srswkt = "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.01745329251994328,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]]";
    // Setup warp options.
    GDALWarpOptions* psWarpOptions = GDALCreateWarpOptions();
    psWarpOptions->dfWarpMemoryLimit = 10 * 1024 * 1024;
    psWarpOptions->hSrcDS = inds;
    psWarpOptions->hDstDS = 0;
    psWarpOptions->nBandCount = 3;
    psWarpOptions->panSrcBands =
        (int*)CPLMalloc(sizeof(int) * psWarpOptions->nBandCount);
    psWarpOptions->panSrcBands[0] = 1;//nessceary and 1 based
    psWarpOptions->panSrcBands[1] = 2;
    psWarpOptions->panSrcBands[2] = 3;
    psWarpOptions->panDstBands =
        (int*)CPLMalloc(sizeof(int) * psWarpOptions->nBandCount);
    psWarpOptions->panDstBands[0] = 1;//nessceary and 1 based
    psWarpOptions->panDstBands[1] = 2;
    psWarpOptions->panDstBands[2] = 3;
    psWarpOptions->pfnProgress = GDALDummyProgress;
    double trans0[6];
    inds->GetGeoTransform(trans0);
    void* transformerArg = GDALCreateGenImgProjTransformer3(
        inds->GetProjectionRef(), trans0,
        srswkt.c_str(), trans);
    //psWarpOptions->pTransformerArg = transformerArg;

    void* hTransformArg2 =
        GDALCreateApproxTransformer(GDALGenImgProjTransform,
            transformerArg,2);
    psWarpOptions->pTransformerArg = hTransformArg2;

    psWarpOptions->pfnTransformer = GDALApproxTransform;// GDALGenImgProjTransform;
    GDALDataset* outds2 = (GDALDataset*)GDALCreateWarpedVRT(inds,
        outxsize, outysize, trans, psWarpOptions);
    outds2->SetProjection(srswkt.c_str());
    //outds2->SetDescription("test-temp3.vrt");
    vector<unsigned char> outdata(256 * 256), outdata1(256 * 256), outdata2(256 * 256);

    int src_x = tilex * 256;
    int src_y = tiley * 256;
    int src_w = 256;
    int src_h = 256;
    if (tilez == 0) src_h = 128;

    outds2->GetRasterBand(1)->RasterIO(GF_Read, src_x, src_y, src_w, src_h, outdata.data(), src_w, src_h, GDT_Byte, 0, 0, 0);
    outds2->GetRasterBand(2)->RasterIO(GF_Read, src_x, src_y, src_w, src_h, outdata1.data(), src_w, src_h, GDT_Byte, 0, 0, 0);
    outds2->GetRasterBand(3)->RasterIO(GF_Read, src_x, src_y, src_w, src_h, outdata2.data(), src_w, src_h, GDT_Byte, 0, 0, 0);
    GDALDestroyWarpOptions(psWarpOptions);
    GDALReleaseDataset(outds2);
    GDALReleaseDataset(inds);
    int64_t ms1 = currentMilliSeconds();

    if (tilez == 0 && tilex == 0 && tiley == 0) {
        string tilename = fname + "_000_II.raw";
        FILE* pf = fopen(tilename.c_str(), "wb");
        fwrite(outdata.data(), 1, outdata.size(), pf);
        fclose(pf);
    }
    

    return ms1 - ms0;

}

void tryOneFileTiles2(string fname) {
    cout << "TryII " << fname << endl;
    int64_t sum_ms = 0;
    int64_t max_ms = 0;
    int64_t min_ms = 9999999;
    int max_x = 0;
    int max_y = 0;
    int max_z = 0;
    int numtiles = 0;
    for (int iz = 0; iz <= 4; ++iz) {
        int numXTiles = pow(2, iz);
        int numYTiles = numXTiles / 2;
        if (numYTiles <= 0) numYTiles = 1;
        for (int iy = 0; iy < numYTiles; ++iy) {
            for (int ix = 0; ix < numXTiles; ++ix) {
                int64_t ms1 = tryATile2(fname, iz, iy, ix);
                if (max_ms < ms1) {
                    max_x = ix;
                    max_y = iy;
                    max_z = iz;
                }
                max_ms = max(max_ms, ms1);
                min_ms = min(min_ms, ms1);
                sum_ms += ms1;
                numtiles++;
                if (numtiles % 500 == 0) cout << numtiles << " ";
                if (iz < 3) {
                    cout << iz << "," << iy << "," << ix << "," << ms1<<"|";
                }
            }
        }
    }
    cout << endl;
    cout << "num " << numtiles << "  "  ;
    cout << "averageII " << sum_ms / numtiles << " ms " ;
    cout << "minII " << min_ms << " ms " ;
    cout << "maxII " << max_ms << " ms " ;
    cout << "maxII zyx " << max_z << "," << max_y << "," << max_x << endl;
    cout << "-------------------" << endl;
}



void tryOneFileTiles(string fname) {
    cout << "Try " << fname << endl; 
    int64_t sum_ms = 0;
    int64_t max_ms = 0;
    int64_t min_ms = 9999999;
    int max_x = 0;
    int max_y = 0; 
    int max_z = 0;
    int numtiles = 0;
    for (int iz = 0; iz <= 4; ++iz) {
        int numXTiles = pow(2, iz);
        int numYTiles = numXTiles / 2;
        if (numYTiles <= 0) numYTiles = 1;
        for (int iy = 0; iy < numYTiles; ++iy) {
            for (int ix = 0; ix < numXTiles; ++ix) {
                int64_t ms1 = tryATile(fname, iz, iy, ix);
                if (max_ms < ms1) {
                    max_x = ix;
                    max_y = iy;
                    max_z = iz;
                }
                max_ms = max(max_ms, ms1);
                min_ms = min(min_ms, ms1);
                sum_ms += ms1;
                numtiles++;
                if (numtiles % 500 == 0) cout << numtiles << " ";
                if (iz < 3) {
                    cout << iz << "," << iy << "," << ix << "," << ms1 << "|";
                }
            }
        }
    }
    cout << endl;
    cout << "num " << numtiles << "  "  ;
    cout << "average " << sum_ms / numtiles << " ms "  ; 
    cout << "min " << min_ms << " ms " ;
    cout << "max " << max_ms << " ms " ;
    cout << "max zyx " << max_z << "," << max_y << "," << max_x << endl;
    cout << "-------------------" << endl; 
}

//返回处理毫秒数，包括从文件系统读取数据时间
//测试3通道数据
int64_t tryATile(string fname,int tilez, int tiley, int tilex)
{
    GDALAllRegister();
    int numXTiles = pow(2, tilez);
    double reso = 360.0 / numXTiles / 256.0;
    double left = -180.0 + tilex * 256 * reso;
    double top = 90.0 - tiley * 256 * reso;

    int64_t ms0 = currentMilliSeconds();
    GDALDataset* inds = (GDALDataset*)GDALOpen(fname.c_str(), GA_ReadOnly);
    int outxsize = 256;
    int outysize = 256;
    double trans[6] = { left , reso , 0 , top , 0, -reso };
    string srswkt = "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.01745329251994328,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]]";
    // Setup warp options.
    GDALWarpOptions* psWarpOptions = GDALCreateWarpOptions();
    psWarpOptions->dfWarpMemoryLimit = 10*1024*1024;
    psWarpOptions->hSrcDS = inds;
    psWarpOptions->hDstDS = 0;
    psWarpOptions->nBandCount = 3;
    psWarpOptions->panSrcBands =
        (int*)CPLMalloc(sizeof(int) * psWarpOptions->nBandCount);
    psWarpOptions->panSrcBands[0] = 1;//nessceary and 1 based
    psWarpOptions->panSrcBands[1] = 2;
    psWarpOptions->panSrcBands[2] = 3;
    psWarpOptions->panDstBands =
        (int*)CPLMalloc(sizeof(int) * psWarpOptions->nBandCount);
    psWarpOptions->panDstBands[0] = 1;//nessceary and 1 based
    psWarpOptions->panDstBands[1] = 2;
    psWarpOptions->panDstBands[2] = 3;
    psWarpOptions->pfnProgress = GDALDummyProgress;
    double trans0[6];
    inds->GetGeoTransform(trans0);
    void* transformerArg = GDALCreateGenImgProjTransformer3(
        inds->GetProjectionRef(), trans0,
        srswkt.c_str(), trans);
    //psWarpOptions->pTransformerArg = transformerArg;

    void* hTransformArg2 =
        GDALCreateApproxTransformer(GDALGenImgProjTransform,
            transformerArg,2);
    psWarpOptions->pTransformerArg = hTransformArg2;

    psWarpOptions->pfnTransformer = GDALApproxTransform;// GDALGenImgProjTransform;
    GDALDataset* outds2 = (GDALDataset*)GDALCreateWarpedVRT(inds,
        outxsize, outysize, trans, psWarpOptions);
    outds2->SetProjection(srswkt.c_str());
    //outds2->SetDescription("test-temp3.vrt");
    vector<unsigned char> outdata(256 * 256),outdata1(256*256),outdata2(256*256);
    outds2->GetRasterBand(1)->RasterIO(GF_Read, 0, 0, 256, 256, outdata.data(), 256, 256, GDT_Byte, 0, 0, 0);
    outds2->GetRasterBand(2)->RasterIO(GF_Read, 0, 0, 256, 256, outdata1.data(), 256, 256, GDT_Byte, 0, 0, 0);
    outds2->GetRasterBand(3)->RasterIO(GF_Read, 0, 0, 256, 256, outdata2.data(), 256, 256, GDT_Byte, 0, 0, 0);
    GDALDestroyWarpOptions(psWarpOptions);
    GDALReleaseDataset(outds2); 
    GDALReleaseDataset(inds);
    int64_t ms1 = currentMilliSeconds();

    if (tilez == 0 && tilex == 0 && tiley == 0) {
        string tilename = fname + "_000.raw";
        FILE* pf = fopen(tilename.c_str(), "wb");
        fwrite(outdata.data(), 1, outdata.size(), pf);
        fclose(pf);
    }

    return ms1 - ms0;
}

int tryWaredVrt()
{
    GDALAllRegister();
    string infilename = "E:/data/MOD13C2/ndvi_rgb_rendered.tif";

    GDALDataset* inds = (GDALDataset*)GDALOpen(infilename.c_str(), GA_ReadOnly);
    
    int outxsize = 256;
    int outysize = 256;
    double trans[6] = { -180, 360.0 / 2 / 256, 0 , 90, 0, -360.0 / 2 / 256 };

    string srswkt = "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.01745329251994328,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]]";
    


    // Setup warp options.
    GDALWarpOptions* psWarpOptions = GDALCreateWarpOptions();
    psWarpOptions->hSrcDS = inds;
    psWarpOptions->hDstDS = 0;
    
    psWarpOptions->nBandCount = 3;// 0 means all bands
    psWarpOptions->panSrcBands =
        (int*)CPLMalloc(sizeof(int) * psWarpOptions->nBandCount);
    psWarpOptions->panSrcBands[0] = 1;//nessceary and 1 based
    psWarpOptions->panSrcBands[1] = 2;
    psWarpOptions->panSrcBands[2] = 3;
    psWarpOptions->panDstBands =
        (int*)CPLMalloc(sizeof(int) * psWarpOptions->nBandCount);
    psWarpOptions->panDstBands[0] = 1;//nessceary and 1 based
    psWarpOptions->panDstBands[1] = 2;
    psWarpOptions->panDstBands[2] = 3;
    psWarpOptions->pfnProgress = GDALTermProgress;

    double trans0[6];
    inds->GetGeoTransform(trans0);
    // Establish reprojection transformer.
    void* transformerArg = GDALCreateGenImgProjTransformer3(
        inds->GetProjectionRef(), trans0,
        srswkt.c_str(), trans);
    psWarpOptions->pTransformerArg = transformerArg;


    psWarpOptions->pfnTransformer = GDALGenImgProjTransform;
    cout << "debug 1" << endl;
    GDALDataset* outds2 = (GDALDataset*)GDALCreateWarpedVRT(inds,
        outxsize, outysize, trans, psWarpOptions);
    cout << "debug 2" << endl;
    outds2->SetProjection(srswkt.c_str());

    cout << "debug 3" << endl;
    outds2->SetDescription("test-temp3.vrt");
    cout << "debug 4" << endl;

    vector<unsigned char> outdata( 256  * 256);
    outds2->GetRasterBand(1)->RasterIO(GF_Read, 0, 0, 256, 256, outdata.data(), 256, 256, GDT_Byte, 0, 0, 0);
    cout << "debug value of 105,55 " << (int)outdata[55 * outxsize + 105] << endl;

    FILE* pf = fopen("temp.raw", "wb");
    fwrite(outdata.data(), 1, outdata.size(), pf);
    fclose(pf);
    cout << "inds cnt 0 " << inds->GetRefCount() << " " << outds2->GetRefCount() << endl;
    // GDALDestroyGenImgProjTransformer(transformerArg); // LineA 问题找到了，根据GDAL 源码 vrtwarped.cpp 
    // GDALDestroyGenImgProjTransformer和GDALReleaseDataset 同时存在报错的问题找到了！
    // 根据GDAL 源码 vrtwarped.cpp ，我发现它构造 warpedVRT dataset之后，使用完后释放资源的时候并没有释放 transformerArg
    // 它代码仅仅释放了 psWarpOptions。因此我怀疑它 transformerArg 传给 warpedVRT 以后就把所属权转给 warpedVRT了，
    // 所以只需要释放一次即可。


    cout << "inds cnt 1 " << inds->GetRefCount() << " " << outds2->GetRefCount() << endl;
    
    cout << "inds cnt 2 " << inds->GetRefCount() << " " << outds2->GetRefCount() << endl;
    
    
    //GDALDestroyGenImgProjTransformer(psWarpOptions->pTransformerArg);
    //psWarpOptions->pTransformerArg = 0;
    GDALDestroyWarpOptions(psWarpOptions);
    cout << "inds cnt 3 " << inds->GetRefCount() << " " << outds2->GetRefCount() << endl;
    GDALReleaseDataset(outds2); // lineB .  这里LineA和LineB要是同时存在就会报异常，单独使用一个就没有问题，why？
    cout << "inds cnt 4 " << inds->GetRefCount() << " " << endl;
    GDALReleaseDataset(inds);


    return 0;
}

int tryWarpedToMem() {
    GDALAllRegister();
    string infilename = "E:/data/MOD13C2/ndvi_rgb_rendered.tif";

    GDALDataset* inds = (GDALDataset*)GDALOpen(infilename.c_str(), GA_ReadOnly);
    int outxsize = 256;
    int outysize = 256;

    GDALDriver* driver = (GDALDriver*)GDALGetDriverByName("MEM");
    GDALDataset* outds = driver->Create("temp.mem", outxsize, outysize, 3, GDT_Byte, 0);
    int xsize2 = outds->GetRasterXSize();
    int ysize2 = outds->GetRasterYSize();
    double trans[6] = { -180, 360.0 / 2 / outxsize, 0 , 90, 0, -360.0 / 2 / outysize };
    outds->SetGeoTransform(trans);
    outds->SetProjection("GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.01745329251994328,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]]");

    // Setup warp options.
    GDALWarpOptions* psWarpOptions = GDALCreateWarpOptions();
    psWarpOptions->hSrcDS = inds;
    psWarpOptions->hDstDS = outds;
    psWarpOptions->nBandCount = 3;// 0 means all bands
    psWarpOptions->panSrcBands =
        (int*)CPLMalloc(sizeof(int) * psWarpOptions->nBandCount);
    psWarpOptions->panSrcBands[0] = 1;//nessceary and 1 based
    psWarpOptions->panSrcBands[1] = 2;
    psWarpOptions->panSrcBands[2] = 3;
    psWarpOptions->panDstBands =
        (int*)CPLMalloc(sizeof(int) * psWarpOptions->nBandCount);
    psWarpOptions->panDstBands[0] = 1;//nessceary and 1 based
    psWarpOptions->panDstBands[1] = 2;
    psWarpOptions->panDstBands[2] = 3;
    psWarpOptions->pfnProgress = GDALTermProgress;

    // Establish reprojection transformer.
    psWarpOptions->pTransformerArg =
        GDALCreateGenImgProjTransformer(inds,
            GDALGetProjectionRef(inds),
            outds,
            GDALGetProjectionRef(outds),
            FALSE, 0.0, 1);
    if (psWarpOptions->pTransformerArg == 0) {
        cout << "psWarpOptions->pTransformerArg is bad." << endl;
        return 12;
    }
    psWarpOptions->pfnTransformer = GDALGenImgProjTransform;

    // Initialize and execute the warp operation.
    GDALWarpOperation oOperation;
    CPLErr err = oOperation.Initialize(psWarpOptions);
    cout << "error " << err << endl;
    oOperation.ChunkAndWarpImage(0, 0,
        GDALGetRasterXSize(outds),
        GDALGetRasterYSize(outds));
    GDALDestroyGenImgProjTransformer(psWarpOptions->pTransformerArg);
    GDALDestroyWarpOptions(psWarpOptions);
    //GDALFlushCache(outds); // no need to flush.

    vector<unsigned char> outdata(outysize * outxsize);
    outds->GetRasterBand(1)->RasterIO(GF_Read, 0, 0, outxsize, outysize, outdata.data(), outxsize, outysize, GDT_Byte, 0, 0, 0);

    cout << "debug value of 105,55 " << (int)outdata[55 * outxsize + 105] << endl;

    GDALClose(outds);
    GDALClose(inds);
    return 0;
}
int tryOutputTifFile() {
    GDALAllRegister();
    string infilename = "E:/data/MOD13C2/ndvi_rgb_rendered.tif";

    GDALDataset* inds = (GDALDataset*)GDALOpen(infilename.c_str(), GA_ReadOnly);


    GDALDriver* driver = (GDALDriver*)GDALGetDriverByName("GTiff");
    GDALDataset* outds = driver->Create("temp.tif", 256, 256, 3, GDT_Byte, 0);
    int xsize2 = outds->GetRasterXSize();
    int ysize2 = outds->GetRasterYSize();
    double trans[6] = { -180, 360.0 / 2 / 256, 0 , 90, 0, -360.0 / 2 / 256 };
    outds->SetGeoTransform(trans);
    outds->SetProjection("GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.01745329251994328,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]]");

    // Setup warp options.
    GDALWarpOptions* psWarpOptions = GDALCreateWarpOptions();
    psWarpOptions->hSrcDS = inds;
    psWarpOptions->hDstDS = outds;
    psWarpOptions->nBandCount = 3;// 0 means all bands
    psWarpOptions->panSrcBands =
        (int*)CPLMalloc(sizeof(int) * psWarpOptions->nBandCount);
    psWarpOptions->panSrcBands[0] = 1;//nessceary and 1 based
    psWarpOptions->panSrcBands[1] = 2;
    psWarpOptions->panSrcBands[2] = 3;
    psWarpOptions->panDstBands =
        (int*)CPLMalloc(sizeof(int) * psWarpOptions->nBandCount);
    psWarpOptions->panDstBands[0] = 1;//nessceary and 1 based
    psWarpOptions->panDstBands[1] = 2;
    psWarpOptions->panDstBands[2] = 3;
    psWarpOptions->pfnProgress = GDALTermProgress;

    // Establish reprojection transformer.
    psWarpOptions->pTransformerArg =
        GDALCreateGenImgProjTransformer(inds,
            GDALGetProjectionRef(inds),
            outds,
            GDALGetProjectionRef(outds),
            FALSE, 0.0, 1);
    if (psWarpOptions->pTransformerArg == 0) {
        cout << "psWarpOptions->pTransformerArg is bad." << endl;
        return 12;
    }
    psWarpOptions->pfnTransformer = GDALGenImgProjTransform;

    // Initialize and execute the warp operation.
    GDALWarpOperation oOperation;
    CPLErr err = oOperation.Initialize(psWarpOptions);
    cout << "error " << err << endl;
    oOperation.ChunkAndWarpImage(0, 0,
        GDALGetRasterXSize(outds),
        GDALGetRasterYSize(outds));
    GDALDestroyGenImgProjTransformer(psWarpOptions->pTransformerArg);
    GDALDestroyWarpOptions(psWarpOptions);

    GDALClose(outds);
    GDALClose(inds);

    return 0;
}